We create this package, mvMonitor, from the foundation laid by Kazor et al (2016). The mvMonitor package is designed to make simulation of multi-state multivariate process monitoring statistics easy and straightforward, as well as streamlining the online process monitoring component.
We begin by simulating an autocorrelated interval-step vector \(t\), where the autocorrelated errors of \(t\) have initial value \(\epsilon_1 \sim \mathcal{U}_{[a,b]}\), where \(a = 0.01\) and \(b = 2\). This \(t\) vector will be sinusoidal with period \(\omega = 7 * 24 * 60\) (signifying a weekly period in minute-level observations) and autocorrelation component \(\phi = 0.75\).
omega <- 60 * 24 * 7
phi <- 0.75
a <- 0.01; b <- 2
mean_t_err <- 0.5 * (a + b) # mean of Uniform(0.01,2)
var_t_err <- (b - a) / 12 # variance of Uniform(0.01,2)
# Create an autocorrelated vector of errors for the random variable t
t_err <- vector(length = omega)
# Initialise the errors so that the overall mean and variance match the sim of
# Kazor et al.
t_err[1] <- rnorm(n = 1,
mean = mean_t_err * (1 - phi),
sd = sqrt(var_t_err * (1 - phi ^ 2)))
for(s in 2:omega){
t_err[s] <- phi * t_err[(s - 1)] +
(1 - phi) * rnorm(n = 1,
mean = mean_t_err * (1 - phi),
sd = sqrt(var_t_err * (1 - phi ^ 2)))
}
# Now for the vector itself
t_star <- vector(length = omega)
t_star[1] <- -cos(2 * pi / omega) + t_err[1]
for(s in 2:omega){
t_star[s] <- -cos(2 * pi * s / omega) +
phi * t_err[(s - 1)] + (1 - phi) * t_err[s]
}
# Now we scale the ts to match the ts from the Unif(0.01,2)
t_star_adj <- ((b - a) * (t_star - min(t_star))) /
(max(t_star) - min(t_star)) + a
Now that we have created the underlying autocorrelated and non-stationary process vector \(t\), we perform some simple confirmatory graphing techniques.
mean(t_star); var(t_star)
## [1] 0.2508151
## [1] 0.5043259
cor(t_star_adj[2:omega], t_star_adj[1:(omega - 1)])
## [1] 0.997075
lmtest::dwtest(t_star_adj[2:omega] ~ t_star_adj[1:(omega - 1)])
##
## Durbin-Watson test
##
## data: t_star_adj[2:omega] ~ t_star_adj[1:(omega - 1)]
## DW = 1.6872, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
We will simulate three features, with each feaure operating under \(k\) different states. We will use \(<x(t), y(t), z(t)>\) instead of the \(<x_1(t), x_2(t), x_3(t)>\) notation used in Kazor et al. This will assist us when we introduce feature states. Our state notation will be \(<x_k(t), y_k(t), z_k(t)>\) for State \(k\).
We create the three features under State 1 (normal operating conditions, or ) as three functions of \(t\): \[\begin{align} \textbf{x}(\textbf{t}) &\equiv \textbf{t} + \boldsymbol\varepsilon_1, \\ \textbf{y}(\textbf{t}) &\equiv \textbf{t} ^ 2 - 3 * \textbf{t} + \boldsymbol\varepsilon_2, \\ \textbf{z}(\textbf{t}) &\equiv -\textbf{t} ^ 3 + 3 * \textbf{t} ^ 2 + \boldsymbol\varepsilon_3, \end{align}\]where \(\varepsilon_i \sim N(0, 0.01)\).
library(scatterplot3d)
library(magrittr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
normal_df <- data.frame(t = t_star_adj,
err1 = rnorm(n = omega, sd = sqrt(0.01)),
err2 = rnorm(n = omega, sd = sqrt(0.01)),
err3 = rnorm(n = omega, sd = sqrt(0.01)))
### State 1 ###
normal_df %<>% mutate(x = t + err1,
y = t ^ 2 - 3 * t + err2,
z = -t ^ 3 + 3 * t ^ 2 + err3)
# plot(normal_df$t)
# plot(normal_df$x)
# plot(normal_df$y)
# plot(normal_df$z)
# # It checks out
orig_state_mat <- normal_df %>% select(x,y,z) %>% as.matrix
Let’s take a look at our 3D scatterplot:
scatterplot3d(x = normal_df$x,
y = normal_df$y,
z = normal_df$x)
Now we can create data for our other two states. These states will be scaled rotations of the current \(<x,y,z>\) set. Thus, we need the rot_and_scale3d() function from the Rotation_Playground.R file. The first state is yaw, pitch, and roll rotated by (0, 90, 30) degrees, and the scales are multiplied by (1, 0.5, 2).
source("inst/Rotation_Playground.R")
### State 2 ###
state2_angles <- list(yaw = 0,
pitch = 90,
roll = 30)
state2_scales <- c(1,0.5,2)
state2_mat <- orig_state_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
We can see the effects of the rotation and scaling.
scatterplot3d(state2_mat)
The second state is yaw, pitch, and roll rotated by (90, 0, -30) degrees, and the scales are multiplied by (0.25, 0.1, 0.75).
### State 3 ###
state3_angles <- list(yaw = 90,
pitch = 0,
roll = -30)
state3_scales <- c(0.25,0.1,0.75)
state3_mat <- orig_state_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
We can see the difference that State 3 creates.
scatterplot3d(state3_mat)
Now, let’s combine all these state observations into one data frame, which we will then subset based on the state indicator.
### Combination and Cleaning ###
# Now concatenate these vectors
normal_df$xState2 <- state2_mat[,1]
normal_df$yState2 <- state2_mat[,2]
normal_df$zState2 <- state2_mat[,3]
normal_df$xState3 <- state3_mat[,1]
normal_df$yState3 <- state3_mat[,2]
normal_df$zState3 <- state3_mat[,3]
# Add a State label column
normal_df %<>%
mutate(state = rep(rep(c(1, 2, 3), each = 60), (24 / 3) * 7))
# Add a date-time column
normal_df %<>%
mutate(dateTime = seq.POSIXt(from = as.POSIXct("2016-11-27 00:00:00 CST"),
by = "min", length.out = 10080))
### Hourly Switching Process ###
state1_df <- normal_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
state2_df <- normal_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
state3_df <- normal_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
normal_switch_df <- bind_rows(state1_df, state2_df, state3_df) %>%
arrange(dateTime)
Let’s take a look at the individual states, and their combined product.
normal_switch_df %>%
filter(state == 1) %>%
select(x, y, z) %>%
scatterplot3d(xlim = c(-2, 3),
ylim = c(-3, 1),
zlim = c(-1, 5))
normal_switch_df %>%
filter(state == 2) %>%
select(x, y, z) %>%
scatterplot3d(xlim = c(-2, 3),
ylim = c(-3, 1),
zlim = c(-1, 5))
normal_switch_df %>%
filter(state == 3) %>%
select(x, y, z) %>%
scatterplot3d(xlim = c(-2, 3),
ylim = c(-3, 1),
zlim = c(-1, 5))
normal_switch_df %>%
select(x, y, z) %>%
scatterplot3d(xlim = c(-2, 3),
ylim = c(-3, 1),
zlim = c(-1, 5))
We now create the faults for individual states and all states. We will introduce these faults 1000 time units after the end of training (to evauluate false positive rates). The first fault is to add 2 to each \(<\textbf{x}(\textbf{t}), \textbf{y}(\textbf{t}), \textbf{z}(\textbf{t})>\).
faultStart <- 8500
### Fault 1A ###
# Shift each <x, y, z> by 2
shift <- 2
fault1A_df <- normal_df %>%
select(dateTime, state, t, x, y, z) %>%
mutate(x = x + shift,
y = y + shift,
z = z + shift)
fault1A_mat <- fault1A_df %>% select(x,y,z) %>% as.matrix
# State 2
fault1A_S2_mat <- fault1A_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
scatterplot3d(fault1A_S2_mat)
# State 3
fault1A_S3_mat <- fault1A_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
scatterplot3d(fault1A_S3_mat)
# Now concatenate these vectors
fault1A_df$xState2 <- fault1A_S2_mat[,1]
fault1A_df$yState2 <- fault1A_S2_mat[,2]
fault1A_df$zState2 <- fault1A_S2_mat[,3]
fault1A_df$xState3 <- fault1A_S3_mat[,1]
fault1A_df$yState3 <- fault1A_S3_mat[,2]
fault1A_df$zState3 <- fault1A_S3_mat[,3]
# Hourly switching process
fault1A_state1_df <- fault1A_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault1A_state2_df <- fault1A_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault1A_state3_df <- fault1A_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault1A_switch_df <- bind_rows(fault1A_state1_df,
fault1A_state2_df,
fault1A_state3_df) %>%
arrange(dateTime)
Now, let’s take a look at our process after introducing fault 1A.
fault1A_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
The second fault (Fault 1B) is to add 2 to feature \(<\textbf{x}(\textbf{t})>\) only.
### Fault 1B ###
# Shift each x by 2
fault1B_df <- normal_df %>%
select(dateTime, state, t, x, y, z) %>%
mutate(x = x + shift)
fault1B_mat <- fault1B_df %>% select(x,y,z) %>% as.matrix
# State 2
fault1B_S2_mat <- fault1B_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
scatterplot3d(fault1B_S2_mat)
# State 3
fault1B_S3_mat <- fault1B_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
scatterplot3d(fault1B_S3_mat)
# Now concatenate these vectors
fault1B_df$xState2 <- fault1B_S2_mat[,1]
fault1B_df$yState2 <- fault1B_S2_mat[,2]
fault1B_df$zState2 <- fault1B_S2_mat[,3]
fault1B_df$xState3 <- fault1B_S3_mat[,1]
fault1B_df$yState3 <- fault1B_S3_mat[,2]
fault1B_df$zState3 <- fault1B_S3_mat[,3]
# Hourly switching process
fault1B_state1_df <- fault1B_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault1B_state2_df <- fault1B_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault1B_state3_df <- fault1B_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault1B_switch_df <- bind_rows(fault1B_state1_df,
fault1B_state2_df,
fault1B_state3_df) %>%
arrange(dateTime)
We now look at the effect of Fault 1B on the scatterplot of the process.
fault1B_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
The third fault we introduce (Fault 2A) is a drift in each feature by \(10 ^ {-3} * (\text{step} - 8500)\).
### Fault 2A ###
# Drift each <x, y, z> by (step - faultStart) / 10 ^ 3
drift_vec <- (1:omega - faultStart) / 10 ^ 3
drift_vec[drift_vec < 0] <- 0
fault2A_df <- normal_df %>%
select(dateTime, state, t, x, y, z) %>%
mutate(x = x + drift_vec,
y = y + drift_vec,
z = z + drift_vec)
fault2A_mat <- fault2A_df %>% select(x,y,z) %>% as.matrix
# State 2
fault2A_S2_mat <- fault2A_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
# scatterplot3d(fault2A_S2_mat)
# State 3
fault2A_S3_mat <- fault2A_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
# scatterplot3d(fault2A_S3_mat)
# Now concatenate these vectors
fault2A_df$xState2 <- fault2A_S2_mat[,1]
fault2A_df$yState2 <- fault2A_S2_mat[,2]
fault2A_df$zState2 <- fault2A_S2_mat[,3]
fault2A_df$xState3 <- fault2A_S3_mat[,1]
fault2A_df$yState3 <- fault2A_S3_mat[,2]
fault2A_df$zState3 <- fault2A_S3_mat[,3]
# Hourly switching process
fault2A_state1_df <- fault2A_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault2A_state2_df <- fault2A_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault2A_state3_df <- fault2A_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault2A_switch_df <- bind_rows(fault2A_state1_df,
fault2A_state2_df,
fault2A_state3_df) %>%
arrange(dateTime)
Now we look at the scatterplot under this drift fault.
fault2A_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
The fourth fault we introduce (Fault 2B) is a drift in features \(\textbf{y}\) and \(\textbf{z}\) by \(10 ^ {-3} * (\text{step} - 8500)\).
### Fault 2B ###
# Drift each <y, z> by (step - faultStart) / 10 ^ 3
fault2B_df <- normal_df %>%
select(dateTime, state, t, x, y, z) %>%
mutate(y = y + drift_vec,
z = z + drift_vec)
fault2B_mat <- fault2B_df %>% select(x,y,z) %>% as.matrix
# State 2
fault2B_S2_mat <- fault2B_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
scatterplot3d(fault2B_S2_mat)
# State 3
fault2B_S3_mat <- fault2B_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
scatterplot3d(fault2B_S3_mat)
# Now concatenate these vectors
fault2B_df$xState2 <- fault2B_S2_mat[,1]
fault2B_df$yState2 <- fault2B_S2_mat[,2]
fault2B_df$zState2 <- fault2B_S2_mat[,3]
fault2B_df$xState3 <- fault2B_S3_mat[,1]
fault2B_df$yState3 <- fault2B_S3_mat[,2]
fault2B_df$zState3 <- fault2B_S3_mat[,3]
# Hourly switching process
fault2B_state1_df <- fault2B_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault2B_state2_df <- fault2B_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault2B_state3_df <- fault2B_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault2B_switch_df <- bind_rows(fault2B_state1_df,
fault2B_state2_df,
fault2B_state3_df) %>%
arrange(dateTime)
We now inspect the effect of Fault 2B.
fault2B_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
Our fifth fault (Fault 3A) is to amplify the underlying \(\textbf{t}\) signal vector as \(\frac{3}{2 * 10080} * (10080 - \text{step}) * \textbf{t}\).
### Fault 3A ###
# Amplify t for each <x, y, z> by 3 * (omega - step) * t / (2 * omega)
amplify_vec <- 2 * (1:omega - faultStart) / (omega - faultStart) + 1
amplify_vec[amplify_vec < 1] <- 1
fault3A_df <- normal_df %>%
select(dateTime, state, t, err1, err2, err3) %>%
mutate(t = amplify_vec * t) %>%
mutate(x = t + err1,
y = t ^ 2 - 3 * t + err2,
z = -t ^ 3 + 3 * t ^ 2 + err3)
fault3A_mat <- fault3A_df %>% select(x,y,z) %>% as.matrix
# State 2
fault3A_S2_mat <- fault3A_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
scatterplot3d(fault3A_S2_mat)
# State 3
fault3A_S3_mat <- fault3A_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
scatterplot3d(fault3A_S3_mat)
# Now concatenate these vectors
fault3A_df$xState2 <- fault3A_S2_mat[,1]
fault3A_df$yState2 <- fault3A_S2_mat[,2]
fault3A_df$zState2 <- fault3A_S2_mat[,3]
fault3A_df$xState3 <- fault3A_S3_mat[,1]
fault3A_df$yState3 <- fault3A_S3_mat[,2]
fault3A_df$zState3 <- fault3A_S3_mat[,3]
# Hourly switching process
fault3A_state1_df <- fault3A_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault3A_state2_df <- fault3A_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault3A_state3_df <- fault3A_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault3A_switch_df <- bind_rows(fault3A_state1_df,
fault3A_state2_df,
fault3A_state3_df) %>%
arrange(dateTime)
Now inspect the effect of this fault
fault3A_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
Our last fault under consideration, Fault 3B, is a dampening fault to \(\textbf{t}\) as \(\log \textbf{t}\). As \(\textbf{t}\) is bdd within \([0.01, 2]\), \(\log\textbf{t}\in[-4.61, 0.7]\).
### Fault 3B ###
# Dampen t for z by log|t|
t_log <- normal_df$t
t_log[faultStart:omega] <- log(t_log[faultStart:omega])
fault3B_df <- normal_df %>%
select(dateTime, state, t, err1, err2, err3) %>%
mutate(t_damp = t_log) %>%
mutate(x = t + err1,
y = t ^ 2 - 3 * t + err2,
z = -t_damp ^ 3 + 3 * t_damp ^ 2 + err3)
fault3B_mat <- fault3B_df %>% select(x,y,z) %>% as.matrix
# State 2
fault3B_S2_mat <- fault3B_mat %*% rot_and_scale3d(rot_angles = state2_angles,
scale_factors = state2_scales)
# scatterplot3d(fault3B_S2_mat)
# State 3
fault3B_S3_mat <- fault3B_mat %*% rot_and_scale3d(rot_angles = state3_angles,
scale_factors = state3_scales)
# scatterplot3d(fault3B_S3_mat)
# Now concatenate these vectors
fault3B_df$xState2 <- fault3B_S2_mat[,1]
fault3B_df$yState2 <- fault3B_S2_mat[,2]
fault3B_df$zState2 <- fault3B_S2_mat[,3]
fault3B_df$xState3 <- fault3B_S3_mat[,1]
fault3B_df$yState3 <- fault3B_S3_mat[,2]
fault3B_df$zState3 <- fault3B_S3_mat[,3]
# Hourly switching process
fault3B_state1_df <- fault3B_df %>%
filter(state == 1) %>%
select(dateTime, state, x, y, z)
fault3B_state2_df <- fault3B_df %>%
filter(state == 2) %>%
select(dateTime, state, x = xState2, y = yState2, z = zState2)
fault3B_state3_df <- fault3B_df %>%
filter(state == 3) %>%
select(dateTime, state, x = xState3, y = yState3, z = zState3)
fault3B_switch_df <- bind_rows(fault3B_state1_df,
fault3B_state2_df,
fault3B_state3_df) %>%
arrange(dateTime)
We inspect the effect of the last fault on our scatterplot.
fault3B_switch_df %>%
select(x, y, z) %>%
scatterplot3d()
The normal_df data frame is without fault. We will remove rows 8500 to 10080, and replace these rows with the same rows from each each of the fault data frames. We then store each of these seven xts matrices in a list.
Implementation notes: turn all the data frames into xts objects for Kazor’s function to work, which we believe follows a depreciated version of zoo, so we need the dateTime column as the rownames of the data frame, not as an index. EDITED NOTE - B. Barnard and I rebuilt Kazor’s code as a functional package called mvMonitoring. The functions in these packages use xts matrices as inputs.
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
normal_switch_xts <- xts(normal_switch_df[,-1],
order.by = normal_switch_df[,1])
# Fault 1A
normal_w_fault1A_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault1A_switch_df[faultStart:omega,])
fault1A_xts <- xts(normal_w_fault1A_df[,-1], order.by = normal_switch_df[,1])
# Fault 1B
normal_w_fault1B_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault1B_switch_df[faultStart:omega,])
fault1B_xts <- xts(normal_w_fault1B_df[,-1], order.by = normal_switch_df[,1])
# Fault 2A
normal_w_fault2A_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault2A_switch_df[faultStart:omega,])
fault2A_xts <- xts(normal_w_fault2A_df[,-1], order.by = normal_switch_df[,1])
# Fault 2B
normal_w_fault2B_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault2B_switch_df[faultStart:omega,])
fault2B_xts <- xts(normal_w_fault2B_df[,-1], order.by = normal_switch_df[,1])
# Fault 3A
normal_w_fault3A_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault3A_switch_df[faultStart:omega,])
fault3A_xts <- xts(normal_w_fault3A_df[,-1], order.by = normal_switch_df[,1])
# Fault 3B
normal_w_fault3B_df <- bind_rows(normal_switch_df[1:(faultStart - 1),],
fault3B_switch_df[faultStart:omega,])
fault3B_xts <- xts(normal_w_fault3B_df[,-1], order.by = normal_switch_df[,1])
faults_ls <- list(normal = normal_switch_xts,
fault1A = fault1A_xts,
fault1B = fault1B_xts,
fault2A = fault2A_xts,
fault2B = fault2B_xts,
fault3A = fault3A_xts,
fault3B = fault3B_xts)
Now that we’ve created and stored these feature, we can use time-series plots to confirm our multi-state vision as well as quickly ascertain the fault behavior.
# Double check that faults are visible and according to plan:
# faults_ls$fault1A %>% select(x) %>% plot.ts()
# faults_ls$fault1A %>% select(y) %>% plot.ts()
# faults_ls$fault1A %>% select(z) %>% plot.ts()
# faults_ls$fault1B %>% select(x) %>% plot.ts()
# faults_ls$fault1B %>% select(y) %>% plot.ts()
# faults_ls$fault1B %>% select(z) %>% plot.ts()
# faults_ls$fault2A %>% select(x) %>% plot.ts()
# faults_ls$fault2A %>% select(y) %>% plot.ts()
# faults_ls$fault2A %>% select(z) %>% plot.ts()
# faults_ls$fault2B %>% select(x) %>% plot.ts()
# faults_ls$fault2B %>% select(y) %>% plot.ts()
# faults_ls$fault2B %>% select(z) %>% plot.ts()
# faults_ls$fault3A %>% select(x) %>% plot.ts()
# faults_ls$fault3A %>% select(y) %>% plot.ts()
# faults_ls$fault3A %>% select(z) %>% plot.ts()
# faults_ls$fault3B %>% select(x) %>% plot.ts()
# faults_ls$fault3B %>% select(y) %>% plot.ts()
# faults_ls$fault3B %>% select(z) %>% plot.ts()
faults_ls$fault1A$x %>% plot.xts()
faults_ls$fault1A$y %>% plot.xts()
faults_ls$fault1A$z %>% plot.xts()
faults_ls$fault1B$x %>% plot.xts()
faults_ls$fault1B$y %>% plot.xts()
faults_ls$fault1B$z %>% plot.xts()
faults_ls$fault2A$x %>% plot.xts()
faults_ls$fault2A$y %>% plot.xts()
faults_ls$fault2A$z %>% plot.xts()
faults_ls$fault2B$x %>% plot.xts()
faults_ls$fault2B$y %>% plot.xts()
faults_ls$fault2B$z %>% plot.xts()
faults_ls$fault3A$x %>% plot.xts()
faults_ls$fault3A$y %>% plot.xts()
faults_ls$fault3A$z %>% plot.xts()
faults_ls$fault3B$x %>% plot.xts()
faults_ls$fault3B$y %>% plot.xts()
faults_ls$fault3B$z %>% plot.xts()
We now throw all this data we’ve just generated into our package. CURRENT ISSUES: remove only alarmed observations, not flagged observations. Set up a replicate structure to repeat this proccess 1000 times.
Sys.time()
## [1] "2017-01-12 15:34:15 MST"
results_ls <- mvMonitoring::mspMonitor(data = faults_ls$normal[,2:4],
labelVector = faults_ls$normal[,1],
trainObs = 2880,
updateFreq = 1440)
Sys.time() # 14 seconds for 720 train, 360 update; 5 seconds for 1440, 720;
## [1] "2017-01-12 15:34:16 MST"
# 1 second for 2880 train (two days), and 1440 update (one day).
falseAlarmRate_SPE <- sum(results_ls$Alarms[,1] == "Alarm") /
nrow(results_ls$Alarms)
falseAlarmRate_T2 <- sum(results_ls$Alarms[,2] == "Alarm") /
nrow(results_ls$Alarms)
faultTime <- zoo::index(faults_ls$normal[faultStart,])
# Why don't we have 10080 - 2880 rows in the Alarm matrix??
str(results_ls$Alarms)
## An 'xts' object on 2016-12-03 00:02:00/2016-12-03 23:59:00 containing:
## Data: chr [1:1434, 1:2] "Normal" "Normal" "Normal" "Normal" ...
## Indexed by objects of class: [POSIXct,POSIXt] TZ:
## xts Attributes:
## NULL